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ABSTRACT 

We study the circularization of tidally disrupted stars on bound orbits around spinning 
supermassive black holes by performing three-dimensional smoothed particle hydrody¬ 
namic simulations with Post-Newtonian corrections. Our simulations reveal that debris 
circularization depends sensitively on the efficiency of radiative cooling. There are two 
stages in debris circularization if radiative cooling is inefficient: first, the stellar debris 
streams self-intersect due to relativistic apsidal precession; shocks at the intersection 
points thermalize orbital energy and the debris forms a geometrically thick, ring-like 
structure around the black hole. The ring rapidly spreads via viscous diffusion, leading 
to the formation of a geometrically thick accretion disc. In contrast, if radiative cool¬ 
ing is efficient, the stellar debris circularizes due to self-intersection shocks and forms 
a geometrically thin ring-like structure. In this case, the dissipated energy can be 
emitted during debris circularization as a precursor to the subsequent tidal disruption 
ffare. The circularization timescale is remarkably long in the radiatively efficient cool¬ 
ing case, and is also sensitive to black hole spin. Specifically, Lense-Thirring torques 
cause dynamically important nodal precession, which significantly delays debris cir¬ 
cularization. On the other hand, nodal precession is too slow to produce observable 
signatures in the radiatively inefficient case. Since the stellar debris is optically thick 
and its photon diffusion time is likely longer than the timescale of shock heating, our 
inefficient cooling scenario is more generally applicable in eccentric tidal disruption 
events (TDEs). However, in parabolic TDEs for Mbh ^ 2 x 10^ Mq, the spin-sensitive 
behavior associated with efficient cooling may be realized. 

Key words: accretion, accretion discs - black hole physics - galactic: nuclei - hy¬ 
drodynamics 


1 INTRODUCTION 


Most galaxies are thought to harbor supermassive black 
holes (SMBHs) with masses from 10^ to 10^ Mq at their 
centers ( [Kormendy Richstone' 1995 Kormendy & Ho 
2Q13| . This is inferred from observing proper motions of 
stars bound to the SMBHs (Schodel et al 2002), measur¬ 


ing stellar velocity dispersions around SMBHs (Magorrian 
et al.||1998 ) or detecting radiation emitted from gas accret¬ 
ing onto the SMBHs (Miyoshi et al. 1995). In the last of 


these signatures, continuous accretion from a gas reservoir 
in active galactic nuclei (AGNs) produces intense radiation 
and powerful outflows and jets. On the other hand, gas ac¬ 
cretion proceeds quiescently, at a significantly lower rate, in 
the centers of inactive galaxies. The gas poor environment 
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surrounding these SMBHs is not accompanied by significant 
emission ( Genzel et al.||2Q03| . 

Tidal disruption events (TDEs) provide a distinctive 
opportunity to probe dormant SMBHs at the centres of such 
inactive galaxies. Most TDEs take place when a star at large 
separation (~ 1 pc) is perturbed onto a parabolic orbit ap¬ 
proaching close enough to the SMBH to be ripped apart 
by the tidal forces, at the radius rt — (Mbh/ui*)^'^^ r* = 
24(MBH/lO®M0)-2/3(m*/M0)-i/3(r*/R0)rs. Here we 
denote SMBH mass with Mbh, stellar mass with m* and 
radius with r*, and the Schwarzschild radius with rs = 
2GMbh/c^, where G and c are Newton’s constant and the 
speed of light, respectively. The subsequent accretion of stel¬ 
lar debris falling back to the SMBH produces a character¬ 
istic flare with a luminosity large enough to exceed the Ed¬ 
dington luminosity for a time scale of weeks to month (Rees 
1988 Evans Kochanek|1989 ). Recent observations of Swift 
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J164449.3+573451 showed that relativistic jets are associ¬ 
ated with some fraction of TDEs ( Bloom et al.||2011 Bur¬ 
rows et al. 2011 Zauderer et al. 2011 Levan et ah 2011). 


Candidates for TDEs have also been observed at X-ray, ul- 


traviolet, and optical wavebands (Komossa V Bade 

1999 

Gezari et a 

.. 1120121 lArcavi et al.||2014| |Holoien et al. 

2014 

Vinko et al. 

2015|), with inferred event rates of 10 ^ per year 

per galaxy ( 

Donley et al.||2002 van Velzen V Earrar||2014 ), 

although the observed light curves and spectra (Gezari et al. 


2012 ) do not always match the simplest theoretical expec¬ 
tations ( Loeb Ulmer|19'^ Lodato Rossi[2011 Strubbe 


Quataert||2011 ). 

Black hole spin is one of two fundamental quantities 
characterizing astrophysical black holes, which inevitably 
acquire spin angular momentum as a result of standard mass 
accretion or chaotic accretion ( King Pringle||2006 ). Mea¬ 
suring black hole spins has proven much more difhcult than 
black hole mass estimation, because the dynamical effects 
of spin occur much closer to the event horizon. Since an ac¬ 
cretion disk can get closer to the black hole when the black 


hole is spinning (Bardeen et al. 1972), a detailed spectral 


analysis of disk X-ray emission can determine the black hole 
spin ( Tanaka et al.||1995 ). Such indirect spin measurements 
have been so far made for about 30 SMBHs ( Miller|[2007 ), 
and recently indicated that the SMBH at the centre of the 
nearby galaxy NGC 1365 has at least 84% of the maximum 
theoretically allowed value ( Risaliti et al.||2013 ). 

SMBH spins are difficult to measure, but are of signifi¬ 
cant astrophysical importance. Spin amplitude and direction 
significantly affect the efficiency for converting rest mass en¬ 
ergy into radiation. While the mass-to-energy conversion ef¬ 
ficiency reaches ^42% for an extreme Kerr black hole in a 
prograde rotation, it is only 4% for the retrograde case 


(Kato et al. 2008), suggesting a wide range of bolometric 


disk luminosities depending on the relative inclination of 
disk and SMBH spin. A SMBH-disk system can also work 
as an engine to convert the black hole’s rotational energy 
into outflows and jets ( Blandford &: Znajek||1977 Koide et 
al.|2002 ). The outflow efficiency depends on spin magnitude 
and direction via a large-scale magnetic flux threading the 
black hole and the disk ( Tchekhovskoy et al.||2012 ). 

It is still theoretically uncertain whether such jets will 
align with the black hole spin axis, with the angular momen¬ 
tum vector of the accretion disk, or with some other aspect 


of the magnetic field geometry (Stone & Loeb 2012a). A 


misaligned accretion disk will undergo differential precession 
due to Lense-Thirring frame dragging. While a geometrically 


thin disk warps by the Bardeen-Petterson effect ( Bardeen 
Petterson[1975 ), a geometrically thick disk can precess as a 
rigid-body rotator, as has been seen in general relativistic 


magneto-hydrodynamic (GRMHD) simulations (Eragile et 
al.|2007 |. Very recent GRMHD simulations have also shown 
that a highly magnetized geometrically thick disk can warp 
due to electromagnetic torques ( McKinney et al.||2Q13 ). 

The present spin of a SMBH records the history of gas 
accretion and mergers with other black holes, and statistical 
samples of SMBH spins encode valuable information on the 


growth history of SMBHs in the universe (Volonteri et al. 


2005 Berti & Volonteri 2008). The SMBHs in most AGN 


are thought to have accreted sufficient gas in their active 
phase to be rotating near the extreme Kerr limit ( Doele- 
m an et aT]|2012 ), although events of randomly oriented gas 


clump accretion might be able to produce black holes that 
rotate much more slowly. In contrast to generally aligned 
prograde accretion in AGN, spinning SMBHs undergoing 
TDEs can rotate either retrograde or prograde with respect 
to inflowing gas, with a full range of possible inclinations 
for the transient accretion disk. TDEs therefore act as nat¬ 
ural laboratories for testing theories about accretion and jet 
launching physics over a full range of prograde and retro¬ 
grade inclination angles. 

Our simulations have focused primarily on tidal disrup- 


tion of stars with eccentric, rather than parabolic (|Evans V 

Kochanek|1989| 

Ayal et al.|2000| |Ramirez & Rosswog|2009| 

Guillochon et al. 

2014|), centre of mass trajectories. Although 


the standard two-body scattering mechanism for generating 
TDEs ( Magorrian V Tremaine|1999 Merritt V Wang|2004 ) 
predicts effectively parabolic trajectories (1 — e* < 
other mechanisms can feed stars to SMBHs at lower eccen¬ 
tricities. Among these non-standard sources of TDEs, the 
most promising are binary SMBHs, a recoil accompanying 
a SMBH merger, and the tidal separation of binary stars. 
Recent numerical simulations have shown that observable 
properties of these “eccentric” TDEs significantly deviate 
from those of standard TDEs; in particular, the rate of 
mass return is substantially increased by being cut off at 
a finite time, rather than continuing indefinitely as a power 


law decay (Hayasaki et al. 2013). Because of their natu¬ 
rally limited dynamic range, simulations of eccentric tidal 
disruption were the first to capture relativistic circulariza¬ 
tion of debris around SMBHs, which is extremely computa¬ 
tionally challenging for the canonical parabolic case. These 
simulations found that circularization is driven by general 
relativistic pericentre shift, which causes shocks to form at 
stream self-intersections. The orbital energy dissipated at 
these self-intersections subsequently circularizes the debris 
into a more compact accretion disk ( Hayasaki et al.||2013 ). 
This is in contrast to past Newtonian simulations of circular¬ 
ization for parabolic orbits around intermediate mass black 
holes (~ 10^ Mq), where purely hydrodynamic effects cir¬ 
cularize tidally stretched debris ( Ramirez V Rosswog|20d9 ); 
although these are expected to be ineffective for SMBH-like 
mass ratios ( Guillochon et al^|2Q14 ). 

The primary motivation for this work is to investi¬ 
gate the effect of SMBH spin on debris circularization. 
We test the hypothesis that nodal precession due to the 
Lense-Thirring effect can delay the onset of stream self¬ 
intersections and strongly retard formation of a luminous 


accretion disk (Gannizzo et al. 


1990 


Kochanek 1994). If 


proven true, TDE circularization delays could be used as 
probes of SMBH spin. Such delays could also decrease the 
average luminosity of many TDEs; if spin-induced circular¬ 
ization delay is common but not universal, it would produce 
a bimodality in TDE optical emission that could explain 
a discrepancy between theoretically predicted and observa- 
tionally inferred TDE rates ( Stone Metzger|2014 ). 

In this paper, we study the circularization of a tidally 
disrupted star on an eccentric orbit around a spinning 
SMBH. In section we describe our numerical approach, 
focusing on the Post-Newtonian corrections we make use 
of. In section we examine the results of our numerical 
simulations in two limiting regimes: one is the radiatively 
efficient cooling case, in which the photon diffusion time is 
much shorter than the energy dissipation timescale. In the 
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opposite scenario, where the radiative cooling is inefficient, 
the debris circularization proceeds in a qualitatively differ¬ 
ent way. In section we examine the effect of black hole 
spin on debris circularization, and the nodal precession of 
the newly formed accretion disk by the Lense-Thirring ef¬ 
fect. Finally, sectionis devoted to summary and discussion 
of our scenario. 


2 METHODS 


We start by describing our numerical methods, with a spe¬ 
cial focus on how to treat relativistic effects in the numerical 
code, and summarize the setup of our physical and numeri¬ 
cal models. First, we describe our procedures for numerically 
modeling the tidal disruption of stars on bound orbits. The 
simulations presented below were performed with a three- 
dimensional (3D) Smoothed Particle Hydrodynamics (SPH) 
code, which is a particle method that divides the fluid into 
a set of particles, and is flexible in setting various initial 
configurations. The code is based on a version originally de¬ 


al. (1995). 


veloped by Benz (1990); Benz et al. (1990); and Bate et 


The SPH equations are composed of a mass conserva¬ 
tion equation, a momentum equation with the SPH standard 
artificial viscosity, and an energy equation. Their details will 
be described later, in section [ tT] These equations, with the 
standard cubic-spline kernel, are integrated using a second- 
order Runge-Kutta-Fehlberg integrator with individual time 
steps for each particle and a variable smoothing length (Bate 


et al. 1995), resulting in enormous computational savings 


when a large range of dynamical timescales are involved. 

The variable smoothing length scheme we used gives 
appropriate spatial resolution in our code, but we ignore the 
term proportional to the gradient of the smoothing length. 
This term is introduced for calculating the gradient of fluid 
properties when the smoothing length is varied in space and 
time, and is important for ensuring energy conservation if 
the gradient of any physical quantities varies over a shorter 


scale than the smoothing length (see Bate 1995 for a review). 
In our simulations, the specific energy is well conserved for 
all the models. This shows that the term plays no crucial 
role in our simulations. 

We have performed 3D SPH simulations self- 
consistently modeling a star from before its entry into the 
tidal sphere up to late times, when the stellar debris has cir¬ 
cularized into a disk. We model general relativistic effects, 
including leading order SMBH spin corrections, by incorpo¬ 
rating Post-Newtonian (PN) forces up to 2PN into the SPH 
code. We have run ten pairs of simulations of tidal disrup¬ 
tion events with different parameters. The common param¬ 
eters through all of simulations are following: m* = ITf©, 
r* = IR©, Mbh = 10®Mq, 7 = 5/3, and a unit of run 
time P* = 27rQ~^ = 27t jGm^ ~ 2.8 hr. The total num¬ 
ber of SPH particles used in each simulation is lOOK, where 
K=1000. We also adopt the standard value of the artificial 
viscosity parameters: <asPH = 1 and ydspH = 2 through all 
the simulations. 

Table 1 summarizes each model. Models 1-4 show the 
eccentric TDEs around non-spinning SMBHs with (e, (3) = 
(0.9,1), (0.8,1), (0.7,1), and (0.7,2). Models 5 and 6 have 
the same simulation parameters as Model 4, except that 



Figure 1. Initial configuration of our simulations. The dashed 
white circle and its central small white dot show the tidal disrup¬ 
tion radius rt and the black hole at the origin, respectively. The 
run time t in units of P* and the number of SPH particles Nqph 
are annotated at the top-left corner and the bottom-right corner, 
respectively. Both the a:-axis and the y-axis are normalized by rt . 
The star is initially located at (0.3rt,0) for Models 1-3, and is 
zoomed into the small square inside the main panel. There, the 
white small arrow indicates the velocity vector of the star. 


the black hole is spinning with spin parameters y = 0.9 for 
Model 5 and x = —0.9 for Model 6. Both models have an 
inclination angle z = 0° between the spin angular momen¬ 
tum and the axis perpendicular to the orbital plane of the 
stellar debris. Models 7 and 8 have the same parameters 
as Model 6 but for i = 90° for Model 7 and i = 45° for 
Model 8, respectively. Model 9 has the same parameters as 
Model 8 but for y = 0.9. Model 10 has the same simula¬ 
tion parameters as Model 1 but for (e, (3) = (0.8, 5), and has 
been performed to compare with Model 2a of |Hayasaki et alT] 
(2013) (see section 3.2). Each of these ten sets of simulation 
parameters has been run twice, with two different equations 
of state (adiabatic and polytropic, discussed in more detail 
in section 1^. 


2.1 Treatment of relativistic effects in SPH 


The formalism of Post-Newtonian (PN) hydrodynamics was 
constructed by Blanchet et al. ( 1990| ) for the approximate 
treatment of relativistic effects in a non-covariant frame¬ 
work. Their formalism is applicable to a moderately rela¬ 
tivistic self-gravitating fluid (with gravitational radiation re¬ 
action, if desired), so long as the PN parameter GMpYtjR(? 
(for a typical spatial scale K) never exceeds 10%. It is not, 
however, simple to implement this formalism into existing 
Newtonian SPH codes. 

Eor a typical TDE with an orbital speed z;, the PN pa¬ 
rameter is estimated to be Oiv^jc^) — 10“^ at the tidal 
disruption radius. The magnitude of the self-gravitating po¬ 
tential and thermal energy of the star can be similarly 
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Figure 2. Orbits of two test particles and one SPH particle 
(from our full disruption simulation with the adiabatic equation 
of state) around a spinning SMBH. The parameters of Model 5 
are adopted for the three particles. Each axis is normalized by the 
tidal disruption radius. The initial positions of the particles are 
located at the center of the small yellow circle. The central white 
point shows the black hole with (x? 0 — (0.9,0°). The solid red 
and dashed white lines show the motion of a SPH particle and a 
test particle in the gravitational potential with Post-Newtonian 
corrections (up to 2PN). The dotted line denotes an orbit of a 
test particle moving in the Kerr metric. 


parameterized to be 0{{v ‘^= 10“® and 
0{cslc^) = 10“^ where Cg is the sound speed (for a stellar 
temperature ~ 10^ K), respectively. These order of magni¬ 
tude estimates show that even the lowest PN order terms 
for stellar self-gravitation and thermal energy can be self- 
consistently neglected, even if up to 2PN precision in the 
black hole’s gravity is desired. Because we only need to mod¬ 
ify the SMBH potential, and can continue to treat hydro¬ 
dynamics and gas self-gravity in a Newtonian fashion, it 
becomes much simpler to implement the PN formalism into 
our SPH code. 

In order to treat approximately relativistic effects such 
as pericentre shift and spin-induced precession, we have 
incorporated acceleration terms corrected by the Post- 
Newtonian approximations into the momentum equation of 
SPH particles. The detailed formulae can be seen in Ap¬ 
pendix A. 


2.2 Initial conditions 

We have performed two-stage simulations: a star is first 
modeled as a polytropic gas sphere in hydrostatic equilib¬ 
rium. The tidal disruption process is then simulated by set¬ 
ting the star in motion through the gravitational field of a 
black hole. In our simulations, the black hole is represented 
by a sink particle with the appropriate gravitational mass 
Mbh. All gas particles that fall within a specified accretion 
radius are accreted by the sink particle. We set the accre¬ 


tion radius of the black hole as the radius of the marginally 


stable orbit for the non-spinning black hole (Bardeen et al. 


1972): Tms — 0.12 rt, in all the models. 


The initial position and velocity of the star is given 
by that of a test particle orbiting around the black hole. 
In the test-particle limit, the specific energy and angular 
momentum with PN corrections are given from equations 
(A7) and (A9) by 


etp 


^ tp 


Ei{rBH — 0,t;BH — 0) 
rrii 

Ji(rm — 0,t;BH — 0) 
rrii 


( 1 ) 

( 2 ) 


where the index i refers to a given SPH particle, and Tbh 
and vbh are the position and velocity vector of the black 
hole particle, respectively. This energy and angular momen¬ 
tum should approximately equal their respective Newtonian 
analogues at a distance far away from the black hole. Given 
an initial position and desired pericenter distance, we nu¬ 
merically solve for an initial velocity vector using the PN 
constants of motion. The initial velocity and position vec¬ 
tor in our simulation models are summarized in Table [U 
The number of initial SPH particles iVspH are lOOK for all 
the models. Figure shows an initial configuration of our 
simulations for Models 1-3. 

Figure shows orbits of two test particles and one SPH 
particle (from our full disruption simulation with the adia¬ 
batic equation of state) around a spinning SMBH. The solid 
red and dashed white lines show the motion of a SPH par¬ 
ticle and a test particle in the gravitational potential with 
Post-Newtonian corrections (up to 2PN). The dotted line de¬ 
notes an orbit of a test particle moving in the Kerr metric. 
The orbit of the test particle in the 2PN potential deviates 
slightly from that of the Kerr metric. It is initially identical 
with the orbit of the SPH particle for the first five orbits, but 
the two diverge afterwards because of hydrodynamic forces 
on the SPH particle. 


2.3 Errors of energy and angular momentum 
conservation 


In order to check convergence of energy and angular momen¬ 
tum conservation with PN corrections, we numerically solve 
the two-body problem with PN corrections in a test particle 
limit by using a fourth-order Runge-Kutta method. While 
the energy and angular momentum are well conserved in the 
case of a circular binary, they oscillate with time in an eccen¬ 
tric binary case (but remain conserved in a time-averaged 
sense). The oscillation amplitude grows with increasing or¬ 
bital eccentricity of the test particle and increasing ratio of 
the tidal disruption radius to pericentre distance of the test 
particle. 

In Table we compare the error levels of energy and 
angular momentum conservation between the test-particle 
simulations and the SPH simulations. Each error level is 
measured by 


Se = 



J - J* 


(3) 


where e and j are the time-averaged and number-averaged 
PN values of specific energy and angular momentum dur¬ 
ing the first ten orbits, while e* and j* are the Newtonian 
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Table 1. Tabulated simulation parameters. The first column shows each simulated model number. The second to seventh columns are 
the penetration factor (3 = rp/rT, the initial orbital eccentricity e*, the initial semi-major axis a*, the radial distance between the black 
hole and the initial position of the star, the specific orbital binding energy of the star e* = —(l/2)/3(l — e*)et where et = GM-Qu/ct — 
1.9 X 10^^ [erg/g], and the black hole spin parameter x (with values between 0 and 1), respectively. The eighth column indicates the 
angle between the black hole spin axis and the axis perpendicular to the orbital plane of the stellar debris. The ninth and tenth columns 
are the periods of the most tightly and loosely bound orbits, respectively (see equations ( |15| > and ( |16| >). The eleventh column shows each 
termination time normalized by P* = 27t— 2.8 hr. The last two columns describe the number of SPH particles at the end 
of simulations for radiatively efficient (A^eff) cind inefficient cooling (A^ineff) cases. 


Model 

d 

e* 

a* [rt] 

ro [rt] 

e* [et] 

X 

i 

^mtb [P*] 

^mlb [P*] 

And [P*] 

Neff 

Nfneff 

1 

1 

0.9 

10 

3 

-0.05 

0 

0° 

11 

44 

100 

99142 

96854 (And = 80) 

2 

1 

0.8 

5 

3 

-0.1 

0 

0° 

4 

13 

100 

99540 

91140 

3 

1 

0.7 

10/3 

3 

-0.15 

0 

0° 

2.2 

6.7 

100 

99980 

85675 

4 

2 

0.7 

5/3 

2.5 

-0.3 

0 

0° 

0.76 

2.3 

40 

99830 

74520 

5 

2 

0.7 

5/3 

2.5 

-0.3 

0.9 

0° 

0.76 

2.3 

40 

99824 

81610 

6 

2 

0.7 

5/3 

2.5 

-0.3 

-0.9 

0° 

0.76 

2.3 

40 

99687 

71148 

7 

2 

0.7 

5/3 

2.5 

-0.3 

-0.9 

90° 

0.76 

2.3 

40 

99805 

72233 

8 

2 

0.7 

5/3 

2.5 

-0.3 

-0.9 

45° 

0.76 

2.3 

40 

99869 

69329 

9 

2 

0.7 

5/3 

2.5 

-0.3 

0.9 

45° 

0.76 

2.3 

40 

99907 

82574 

10 

5 

0.8 

1 

1.8 

-0.5 

0 

0 

0.35 

1.03 

10 

99632 

- 


Table 2. Initial conditions and errors for our simulations. The first column shows each simulated model. The sec ond and th ird columns 
denote the initial position and velocity vector for each model. The normalization of the velocity is given by vt = yj GMbh/p- The fourth 
and fifth column show the energy conservation error and angular momentum conservation error of a test particle, respectively. The last 
two columns describe the energy conservation error and angular momentum conservation error of a SPH particle, respectively. 


Model 

ro [rt] 

ro [rt] 

Setp [%] 

'5itp [%] 

(iesPH [%] 

(^jsPH [%] 

1 

(0.0,-3.0, 0.0) 

(0.447,0.589,0.0) 

1.0 

0.08 

0.054 

0.041 

2 

(0.0,-3.0, 0.0) 

(0.436,0.512,0.0) 

1.0 

0.07 

0.041 

0.026 

3 

(0.0,-3.0, 0.0) 

(0.506,0.474,0.0) 

0.27 

0.02 

0.16 

0.05 

4 

(0.0,-2.5, 0.0) 

(0.359,0.25,0.0) 

4.0 

0.45 

0.25 

0.02 

5 

(0.0,-2.5, 0.0) 

(0.359,0.251,0.0) 

3.5 

0.35 

0.12 

0.31 

6 

(0.0,-2.5, 0.0) 

(0.359,0.248,0.0) 

4.8 

2.8 

0.38 

0.23 

7 

(0.0,-2.5, 0.0) 

(0.359,0.249,0.0) 

4.2 

0.45 

0.25 

0.4 

8 

(0.0,-2.5, 0.0) 

(0.359,0.251,0.0) 

4.5 

2.0 

0.19 

0.19 

9 

(0.0,-2.5, 0.0) 

(0.359,0.249,0.0) 

3.5 

0.18 

0.31 

0.6 

10 

(0.556, -1.71,0) 

(0.377,0.1225,0) 

10.6 

1.74 

5.1 

22.6 


specific energy and angular momentum of the initially ap¬ 
proaching star, respectively. Except for Model 10, the energy 
and angular momentum is well conserved at an error level 
of 2%. 


3 STELLAR DEBRIS CIRCULARIZATION 


Recent numerical simulations have shown that the peri- 
centre shift of the stellar debris plays an essential role in 
quickly forming an accretion disk around a non-spinning 
SMBH, because it leads to debris orbit self-intersections, 
which dissipate energy in shocks and cause rapid circulariza¬ 
tion ( Hayasaki et al.||2Q13 ). This work also showed that the 
angular momentum of the stellar debris is conserved during 
circularization. This angular momentum conservation allow 
us to estimate the circularization radius of the stellar debris, 
which is given by 


2x 1 + e* 

rc = a*(l - e*) = —^—rt, 


/3 


(4) 


where a*, e*, and /3 are the semi-major axis and orbital 
eccentricity of the initially approaching star, and the ratio of 
tidal disruption radius n to pericenter distance rp = a*(l — 
e*), respectively. The specific binding energy of the stellar 
debris measured at the circularization radius can then be 
written as 


Cc 


1 /3 

2rTZ^'’ 


(5) 


where et = G^Mbh/p is a characteristic specific energy of 
the tidal disruption radius. On the other hand, the specific 
orbital energy of the initially approaching star is: 

e* =(6) 

Note that the IPN, 1.5PN, and 2PN order terms are propor¬ 
tional to (GMjrt)!(? ~ 2.1%, (GM!(? ~ 0.31%, and 
(GMjrtY ^ 0.045%, respectively, at the tidal disrup¬ 
tion radius. This shows that equations are typically 

corrected at the ~ 2.5% level by our PN approaches. 

The difference between ?Ti*e* and gives the maxi- 
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Figure 3. Evolution of the specific angular momentum and energy for Models 1-9 in the radiatively efficient cooling cases. These constants 
are averaged out per SPH particle. In panels (al)-(a3), the black, blue, and red solid lines denote the specific angular momentum 
normalized by jt = VGMbh^- The corresponding dashed lines show the specific angular momentum of a test particle moving in 
a Newtonian potential, j* = —ef). In panels (bl)-(b3), the black, blue, and red solid lines represent the specific binding 

energy normalized by et = GM-Qa/'f't- The corresponding dashed lines show the Newtonian specific binding energy of a test particle, 
e* = —(1/2)^(1 — e*)et. The dotted line shows the Newtonian specific binding energy measured at the tidal disruption radius. The run 
time t is in units of P* = 21^/ GM ~ 2.8 hr. 


mum amount of binding energy potentially dissipated during 
debris circularization: 

c- _ I I _ m* 1 n V/ 10^2 r ] 

oemax — ?n*|€* Cel — ^ ~ ; T£t — 1.9 X 10 [ergj 


2(1 + e*) 

4/3 / \ -1 


2/3 


X .(7) 


(l + e)\^Mo/ \Rq) \\(f>MQ) 

It is crucial to consider where the dissipated energy goes 
during debris circularization. The photon diffusion timescale 
of the stellar debris is given by ( Mihalas h Mihalas[I984 ) 


( Mbh a ^ ^ 

I 10®M© ) \rt) \ n ) 


where Kes = 0.4 [cm^ g is the opacity for electron scatter¬ 
ing, and 

Eo = ++ 6.5 X 10® [gem-"] ^ 

2nrAr ■' \Rq J \Mq ^ 

( Mbh \ / r/Ar\. 

^ (lOSMo) (rtJ (rt) ^ 


^diff = 


X 


H 





/ m* \ ^ f Mbh 
\M^J VlO^ 



(8) 


where k, E, and r are the scale height, opacity, surface 
density, and optical depth of the stellar debris, respectively. 
The optical depth is approximately estimated to be 


T 


Hi pH 


kJ: - 2.6 X 10® 



is the fiducial surface density, where r and Ar are the ra¬ 
dial size and width of the debris ring, respectively. We note 
that the stellar debris is clearly optically thick. If the pho¬ 
ton diffusion timescale is longer than the energy dissipation 
timescale (i.e. the debris circularization timescale), then ra¬ 
diative cooling is inefficient and dynamically unimportant. 
Otherwise, the radiative cooling can be efficient; the dynam¬ 
ical effect of this will be to reduce the thickness of the debris 
streams. 

^ It is clear from equation Q that in the eccentric TDEs 
we simulate, tdiff will always be very long compared to or¬ 
bital timescales, and cooling will generally be radiatively 
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Figure 4. Evolution of the specific angular momentum and energy for Models 1-9 in the radiatively inefficient cooling cases. The figure 
formats are the same as Figure 


inefficient. However, this is not necessarily true for the 
parabolic (e* ~ 1) TDEs which dominate the event rate. 
If we approximate the frozen-in specific energy spread of 
e* 1 tidal debris as ( Stone et al.|2Q13 ) 

GMbu 


Aet = 


n n 


( 11 ) 


then we can calculate the fallback time with Kepler’s third 
law for the most tightly bound debris to be 


^fb — 


X 


p, = 3.5x10« 

2v 2 V m* / 



Mbh 

WMf 


1/2 


( 12 ) 


We can estimate the regime where cooling is important by 
requiring tdiff < ^fb; this is conservative because circular¬ 
ization likely takes several fallback times to complete. If we 
assume roughly cylindrical debris streams, with H = Ar, 
then factors of Ar cancel and we are left with a simple con¬ 
dition on the maximum extent of the debris stream, r, which 
we hereafter identify as debris apocenter r^. Specifically, 


— > 1.7x10^ 

n 


f ^ \ f ^BH 
wJ V 106Mo 


- 5/6 


Rg) 


\ / mC\ 


7/3 


Mo 7 


This illustrates why eccentric TDEs should generally be 
in the radiatively inefficient limit, but if we substitute in 
the apocenter of the most tightly bound debris streams for 


parabolic TDEs, we obtain the following condition for ra¬ 
diatively efficient cooling: 


-^BH ^ l.GxlO^Afo 





.(14) 


Under our simplifying assumptions, radiative cooling will 
be, in general, likely to play some dynamical role for de¬ 
bris streams of parabolic TDEs. The radiative efficiency of 
these streams could be reduced if they cool to the point 
where bound-free absorption dominates electron scattering 
as a source of opacity (Kochanek 1994), but it could also 
be increased if magnetically driven turbulent advection of 
photons enhances cooling rates ( Jiang et al.|2014 ). 

Given the many uncertainties in this discussion, and the 
possible applicability of both regimes to parabolic TDEs, we 
present two extreme cases: one involves a set of radiatively 
efficient cooling simulations, where the entropy remains con¬ 
stant through the simulation, in the polytropic equation of 
state with 7 = 5/3. The other involves the radiatively inef¬ 
ficient cooling simulations, where the entropy is locally in¬ 
creased with adiabatic equation of state but the total energy 
is conserved. The detailed parameters of ten (9+1) simula- 
•(13)tion models are shown in Table [T] Note that Model 10 is 
done for the purpose of comparison to the fiducial simu¬ 
lation model of Hayasaki et al. (2013), where the pseudo- 
Newtonian potential was adopted. Further details are de¬ 
scribed in section 1321 
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3.1 Radiatively efficient cooling cases 


First, we describe the results of our radiatively efficient cool¬ 
ing simulations, which serve as one extreme of possible ra¬ 
diative cooling. Figure shows the evolution of the spe¬ 
cific angular momentum and specific energy in Models 1- 
9, which are averaged per SPH particle. Panels (bl)-(b3) 
show the evolution of the specific binding energy in Mod¬ 
els 1-9. The thermal energy is estimated to be the order of 
10“^ et for all the models, and therefore it is negligibly small 
compared with the magnitude of the orbital binding energy. 
Panels (al)-(a3) show the evolution of specific angular mo¬ 
mentum per SPH particle. The small-amplitude oscillations 
seen there are due to the PN terms added in the SPH equa¬ 
tion of motion. Equivalent errors are also seen in energy and 
angular momentum conservation for a test particle on an 
eccentric orbit moving under a gravitational potential with 
PN corrections. We compared these errors between the SPH 
simulations and the test particle integrations. The detailed 
results can be seen in columns 4-7 of Table 2 in Section [231 
Note that the angular momentum of the SPH particles in all 
the models is conserved at a less than 2% error level through 
the simulations. 

Figures [5][7| show a sequence of snapshots of the sur¬ 
face density of stellar debris, which is projected on the x-y 
plane in a logarithmic scale, covering two orders of mag¬ 
nitude, for Models 1-3. Each figure progresses from panel 
(a) to panel (d) in chronological order. The central small 
point, dashed circle and white small arrows show the black 
hole, tidal disruption radius, and velocity field of the stellar 
debris, respectively. The run time is noted at the top-right 
corner in units of P*, while the number of SPH particles is 
indicated at the bottom-right corner. 

The stellar debris moves around the black hole for sev¬ 
eral orbits. Over time, the debris stretches due to the spread 
in its constituent orbital energies, and the debris head in¬ 
teracts with the tail near apocentre, leading to significant 
energy dissipation in shocks. The binding energy of the stel¬ 
lar debris is substantially reduced by a sequence of orbit 
crossings, causing the debris to circularize. Erom Eiguresj^ 
and the stellar debris clearly circularizes in Models 1 and 
2. However, circularization has proceeded much less rapidly 
in Model 3, as we see from Eigurej^ 

Panel (bl) shows the evolution of specific binding en¬ 
ergy in Models 1-3. Since the specific binding energy has not 
reduced from e* to Cc even at the end of these runs, the circu¬ 
larization process has not yet completed. Adopting a simple 
extrapolation from e* to Cc, the circularization timescales 
can be estimated to be ~ 120P* in Model 1, ~ 180P* in 
Model 2, and ~ 2500P* in Model 3, respectively. These 
extrapolated timescales indicate strongly varying per-orbit 
efficiencies of shock dissipation: if this efficiency were con¬ 
stant, then tc oc . In fact, the circularization timescale 
behaves in an inverse manner. The prominently long cir¬ 
cularization timescale of Model 3 shows this counterintu¬ 
itive behavior. This declining dissipation efficiency at fixed 
/3 and decreasing e* is likely because the relative velocity 
between the debris head and the debris tail at their self¬ 
intersection decreases as we go from Model 1 to Model 3. 
This is a general feature of eccentric TDEs at fixed /3: low 
eccentricity produces self-intersections closer to apocentre, 
with lower relative velocities. This is also confirmed by the 



Figure 5. A sequence of snapshots of the tidal disruption process 
in Model 1 (a* = 10rt, e* = 0.9, ^ = 1, x = 0.0, and i = 0°) in 
the radiatively efficient case. Panel (a) to (d) are shown in chrono¬ 
logical order. All the panels show a surface density projected on 
the x-y plane for 0 ^ t ^ 100. The color bar shows the magnitude 
of the density in a logarithmic scale, where Eq = 6.5 X 10® gcm“^ 
is the fiducial surface density (see equation |10[ ). The black hole 
is set at the origin. The run time t is in units of P* and is anno¬ 
tated at the top-right corner. The number of SPH particles are 
indicated at the bottom-right corner. The white small arrows and 
the dashed circle indicate the velocity field of the stellar debris 
and the tidal disruption radius, respectively. 


difference between the orbital periods of most tightly and 
loosely bound gas, which shrinks as the orbital eccentricity 
decreases. These periods are given by ( Hayasaki et al.|2013 ) 

1 3/2 


tmtb — 


2V2 


^(1 -e.) 




(15) 


Imlb — 


2V2 


/?(l-e.) 1/3 


- 3/2 


(16) 


where g* = ?ti*/Mbh. Panel (bl) shows that the energy dis¬ 
sipation rate increases with orbital eccentricity for the case 
of /3 = 1. This implies that stellar debris should efficiently 
circularize in most standard, parabolic TDEs (/3 = 1 and 
e* = 1) around non-spinning SMBHs. 


3.2 Radiatively inefficient cooling cases 

Next, we describe the results of our radiatively inefficient 
cooling simulations for Models 1-3. Eigures |8|10| show a se¬ 
quence of snapshots of debris surface density for Models 1-3, 
with the same figure formats as in Eigure.j^ but for the ra¬ 
diatively inefficient cooling simulations. 

In Eigurej^ the stellar debris moves away from the black 
hole following tidal disruption, as shown in panel (a), and is 
then stretched during pericentre return in panel (b). Near 
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Figure 6. A sequence of snapshots of the tidal disruption process 
in Model 2 (a* = 5rt, e* = 0.8, (3 = 1, x = ^-nd i = 0°) in 
the radiatively efficient case. The figure formats are the same as 
Figure 


Figure 8. A sequence of snapshots of the tidal disruption process 
in Model 1 (a* = 10rt, e* = 0.9, ^ = 1, x = ^^^d i = 0°) in 

the radiatively inefficient case. The figure formats are the same 
as Figure but for 0 ^ t ^ 80. 



Figure 7. A sequence of snapshots of the tidal disruption process 
in Model 3 (a* = lOr^/s, e* = 0.7, ^ = 1, x = ^-nd i = 0°) in 
the radiatively efficient case. The figure formats are the same as 
Figure 


Figure 9. A sequence of snapshots of the tidal disruption process 
in Model 2 (a* = 5rt, e* = 0.8, /3 = 1, x = O-O? cind i = 0°) in 
the radiatively inefficient limit. The figure formats are the same 
as Figure but for 0 ^ t ^ 80. 


apocentre, the leading “head” of the debris significantly in¬ 
tersects with the trailing “tail,” as can be seen in panel (c). 
After several tens of orbits, the debris expands significantly, 
as the thermal energy increases from shock heating. This 
can be seen in panel (d). 

There is a clear difference between Model 1 and Models 


2-3. While no accretion disk forms by the end of simulation 
in Model 1, an accretion disk clearly forms and viscously 
evolves in Models 2 and 3. This viscous evolution can be un¬ 
derstood as follows: the viscous timescale for an a-viscosity 
disk is given by 
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Figure 10. A sequence of snapshots of the tidal disruption pro¬ 
cess in Model 3 (a* = 10/3 n, e* = 0.7, l3 = 1, x = O-O? ^^nd 
z = 0°) in the radiatively inefficient limit. The figure formats are 
the same as Figure but for 0 ^ t ^ 80. 


r“ 5 

ly TT \ass. 


^vis ^ ^ 


3/2 /i7\-3/2 


/ 0.1 \ / r y/2 /H\ 
Vttss J \rtJ \ r J 


P*. 


(17) 


The viscous timescale can be comparable to the debris or¬ 
bital period because of the enhanced pressure and resultant 
geometrically thick structure {H/r ~ 1) in the radiatively 
inefficient regime, meaning that rapid viscous redistribution 
of angular momentum and energy is possible. However, our 
simulations become significantly less reliable at late times, 
after formation of the accretion disk, due to the lack of ra¬ 
diation pressure and magneto-hydrodynamics in our code. 
We present these late-time results for completeness, but em¬ 
phasize that our simulations are primarily designed only to 
simulate the circularization process itself. 

Figureshows evolution of the specific angular momen¬ 
tum and specific energy in Models 1-9, which are averaged 
per SPH particle. Panels (dl)-(d3) show the evolution of the 
specific energy in Models 1-9. The specific energy increases 
with time, in stark contrast with the radiatively efficient 
case. Preferential accretion of highly bound SPH particles in¬ 
creases, over time, the mean specific energy of non-accreted 
SPH particles. Panels (cl)-(c3) show the evolution of specific 
angular momentum averaged over all remaining SPH parti¬ 
cles. This also increases with time because of preferential 
accretion of the lowest angular momentum SPH particles. 

As a check on our radiative efficiency assumptions, we 
compare the photon dihusion timescale with the timescale 
for energy dissipation by shock heating. We estimate this lat¬ 
ter timescale using the diherence between the orbital periods 
of the most tightly and loosely bound gas, which are given by 


equations (15) and (16), respectively. These timescales are 


shown for each model in Table ^ Model 1 has the longest 
energy input timescale among all models, on the order of 
10^ s. From Figure 1]_ the surface density of the stellar de¬ 
bris and its scale height are estimated to be E/Eq ~ 10“°'® 


and H/vt ~ 1 for Model 1. Substituting into equation ([^, 
the photon dihusion timescale is ~ 10^ s, clearly longer than 
the shook heating timescale, ft is therefore clear that eccen¬ 
tric TDEs operate in the radiatively inefficient regime, al¬ 
though as we have argued above, the parabolic case is more 
ambiguous. 


3.3 Comparison to pseudo-Newtonian potential 
simulation 


In this section, we compare our SPH simulation with 
2PN corrections to our SPH simulation with a pseudo- 
Newtonian potential. The PN simulation of Model 10, whose 
parameters can be seen in Table 1, was performed for the 
purpose of comparing with the pseudo-Newtonian simula¬ 
tion model. In Model 10, we used the same initial condition 


as that of Model 2a in our previous paper (Hayasaki et al, 


2013). The initial position and velocity are set by 


ro = (ro cos 0, ro sin 0, 0), 

Vo = (r(ro) COS0O - ^o0(0o) sin^o, r(ro) sin^o 

+ ro0(0o) COS0O, 0), 


where ro = aHc(l —6=,= ) and 0o = —0.27r are adopted. Here, the 
radial velocity r and angular velocity r and 0 are given by 
energy conservation and angular momentum conservation as 


r 

0 


2(€p3eudo - f7(r)) - 
^pseudo 


where Cpseudo and /pseudo are the specific energy and the 
specific angular momentum for bound orbits, the pseudo- 
Newtonian potential (Wegg 2012), respectively. They are 
written by 


t/(r) = - 


GMb 


1 — Cl rs 

1 — C 2 (rs /2r) 2r 


Cpsuedo — 


(rp/ra)^[/(rp) - [/(ra 


(rp/ra)2 - 1 

/pseudo = v^2r^(e - U{rp)) = - U{ra)), 


(18) 


where rp = a*(1 — 6=,= ) and ra = a=,c(l + e=,c) are the pericenter 
distance and the apocenter distance, respectively, and we 
adopt that ci = (—4/3) (2 + \/6), C 2 = (4\/6 — 9), and cs = 
(—4/3)(2v^ —3). The initial position and velocity vector for 
Model 10 are seen in Table |2l 

Figure shows the evolution of specific binding en¬ 
ergy and angular momentum. In panel (a), the blue solid, 
red solid, and black dashed lines show the specific angu¬ 
lar momenta of Model 10, our pseudo-Newtonian simulation 
model, and /pseudo respectively. From the panel, we note 
that the angular momentum is conserved for the pseudo- 
Newtonian case, while it is shifted at the ~ 20% level from 
the Newtonian specific angular momentum and at the 3% 
level even from the pseudo-Newtonian case. In panel (b), the 
blue solid, red solid, and black dashed lines show the spe¬ 
cific binding energies of Model 10, our pseudo-Newtonian 
simulation model, and Cpseudo respectively. From the panel, 
the specific binding energy first agrees well with that of test 
particle, but substantially reduces due to debris circular¬ 
ization, and eventually saturates at ~ 7P* in the pseudo- 
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Figure 11. Density map for Model 1 (a* = 10rt, e^c = 0.9, /3 = 1, x = 0.0, and i = 0°) in the radiatively inefficient cooling limit. Each 
panel shows surface densities projected on x-y plane (left panel) and on x-z plane (right panel) at t = 80. The figure formats are the 
same as Figure 


Newtonian case. This is because there is less shock energy 
dissipation required for accretion disk formation (see Figure 
12 of Hayasaki et al.||2013 ). On the other hand, the specific 
binding energy of Model 10 is slightly shifted from that of 
test particle: Cpseudo = e* = — 0.5et, with some small oscilla¬ 
tion. 


The errors of Model 10 in those two panels are at¬ 
tributed to the intrinsic error of PN corrections. In this 
comparison, the error of angular momentum is substantially 
larger than that of binding energy. This trend is different 
from our other simulation results and test particle simula¬ 
tions (see Table 2). This is probably attributed to the higher 
/3 value of Model 10. 


Debris circularization is significantly more efficient in 
the pseudo-Newtonian simulation. This is because the pre¬ 
cession rate is larger than that of Model 10. Since the peri- 
centre velocity of the debris in these models is a significant 
fraction of the speed of light, the PN approximation should 
break down. Thus, our simple PN approaches to relativistic 
effects should not be applicable for such a high 13 simulation 
as /3 = 5. 


Bonnerot et al. (2015) have recently performed similar 


SPH simulations for accretion disk formation around a non¬ 
spinning SMBH. Some of their simulations have used the 


same simulation parameters and setup as (Hayasaki et al. 


2013) but for a recently derived pseudo-Newtonian poten¬ 


tial, isothermal equation of state, and 500K SPH particles. 
We have compared the black line of their Figure 3 (Model 
RI5e.8) with two lines of panel (b) of Figure 12. The specific 
energy evolution of their model, while similar to our older 
pseudo-Newtonian simulation, differs significantly from our 
PN Model 10. Specifically, they see more rapid circulariza¬ 
tion. This could be due to the somewhat different equa¬ 
tions of state employed (isothermal vs polytropic) or our 


lower particle resolution, but is most likely due to differ¬ 
ences in the gravitational potential employed. While both 
potentials accurately account for apsidal precession around 
Schwarzschild SMBHs in /3 = 1 and 13 — 2 events, our PN 
approach begins to break down for Model lO’s /3 = 5, under¬ 
estimating the true precession rate and therefore generating 
an artificially weak (and delayed) stream self-intersection. 
We hope to perform more detailed comparisons in future 
work. 

For ^ — 1 and e = 0.8, our radiatively efficient simu¬ 
lation for Model 2 is the same as their Model Rile.8, but 
for the pseudo-Newtonian potential, isothermal equation of 
state, and 500K SPH particles. We have compared the red 
dashed line of their Figure 10 (Model Rile.8) with the blue 
line of panel (bl) of our Figure 3 at t ~ 90, which corre¬ 
sponds to the end of their simulation. The onset of the en¬ 
ergy dissipation in their model is similar to ours, but their 
dissipated energy estimated at t ~ 90 is ~ 25% smaller than 
our dissipated energy. This could mainly originate from the 
low resolution of the simulation. We will discuss the numer¬ 
ical convergence problem in the next section. 

3.4 Numerical convergence 

In order to test the numerical convergence of our simula¬ 
tions, we have repeated Model 4 with 200K and 500K SPH 
particles (as opposed to our standard lOOK SPH particles). 
Figure shows the evolution of specific angular momen¬ 
tum and energy for Model 4 in both the radiatively efficient 
and inefficient cooling regimes. In panels (a)-(d), the red 
and black lines are the 200K and lOOK cases, respectively. 
Note that the blue solid line of panel (b) shows normalized 
specific binding energy until t ^ 3S for the 500K particles’ 
simulation. The specific angular momentum of the higher 
resolution case is better-conserved, but both simulations are 
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Figure 12. Evolution of the specific angular momentum and specific binding energy in Model 10 and our pseudo-Newtonian simulation 
model. These are averaged out per SPH particle. The specific binding energy and specific angular momentum are normalized by et = 
GM-Qu/'f't and jt = \/GMbh ^5 respectively. In panels (a), the blue and red solid lines denote the specific angular momentum of Model 
10 and pseudo-Newtonian model, respectively. The dashed and dotted lines denote /pseudo and the specific angular momentum of a test 
particle moving under the Newtonian potential, j* = y^a*(l — e*), respectively. In panel (b), the blue and red solid lines represent the 
specific binding energy of Model 10 and pseudo-Newtonian model, respectively. The dashed lines show Cpseudoj which approximately 
equals to the Newtonian specific binding energy of a test particle, e* = —(1/2)^(1 — eHc)et. 


well-converged at a less than 0.2% level for the radiatively 
efficient regime. 


Because the angular momentum is conserved, the spe¬ 
cific binding energy at the circularization radius is estimated 
to be Cc ~ 0.588 for Model 4 from equation The 200K 
simulation achieves complete circularization by the end of 
the simulation, but the lOOK and 500K simulations do not. 
There is a ~ 17% diherence of the final energy distributions 
between lOOK and 200K particles’ simulations, and a ~ 14% 
diherence between the lOOK and 500K simulations. This in¬ 
dicates that the late stages of debris circularization are not 
numerically converged in our radiatively efficient simulations 
with A^sph = lOOK. The 200K case is slightly more dissi¬ 
pative, by ~ 3%, compared with 500K case, because the 
lower resolution makes the artificial viscosity larger. Thus, 
the 500K case is better converged than 200K case. 


An equivalent convergence test on Model 4 with a radia¬ 
tively inefficient (adiabatic) equation of state finds a smaller, 
~ 10 %, diherence between the final energy distributions of 
the lOOK and 200K cases (see panel (d) of Figure 13). Im¬ 
portantly, the radiatively inefficient runs appear much more 
closely converged through the end of the circularization pro¬ 
cess, and they only diverge from each other as the thick torus 


begins to viscously accrete (which is not the main focus of 
this work). Although our radiatively inefficient simulations 
appear to have converged in their description of debris circu¬ 
larization, our radiatively efficient runs are underestimating 
the true efficiency of debris circularization due to resolution 
ehects. This is a topic we will examine more clearly in future 
work. 


4 STELLAR TIDAL DISRUPTION BY A 
SPINNING BLACK HOLE 

Finally, we investigate the tidal disruption of a star by a 
spinning SMBH. First, we see how black hole spin ahects 
coplanar debris circularization through simulations of Mod¬ 
els 4-6. Second, we examine Lense-Thirring precession of 
stellar debris during the tidal disruption process, in Mod¬ 
els 7-9. We have used Model 4 specifically to check our re¬ 
sults for numerical convergence, by running two higher res¬ 
olution versions (one for each equation of state) with 200K 
SPH particles. We find that Model 4 is well-converged for 
the radiatively inefficient regime during the period of debris 
circularization, but in the radiatively efficient regime, our 
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Figure 13. Evolution of the specific angular momentum and energy for Model 4 with A^sph = 200K (as opposed to our standard lOOK 
SPH particles). These are averaged out per SPH particle. Panels (a) and (b) are for the radiatively efficient cooling simulation, whereas 
panels (c) and (d) are for the radiatively inefficient cooling simulation. In panels (a) and (c), the black and red solid lines denote the 
specific angular momentum normalized by jt = for lOOK and 200K SPH particles’ simulations, respectively. In panels (b) and 

(c), the black and red solid lines represent the specific binding energy normalized by et = GM-Qn/rt for lOOK and 200K SPH particles’ 
simulations, respectively. In panel (b), the blue solid line shows the evolution of the normalized specific energy of the A^sph = 500K 
simulation for Model 4. All other figure formats are the same as Figure 


lower-resolution fiducial runs may be underestimating the 
efficiency of shock dissipation (see panel (b) of Figure p^. 


4.1 Effect of black hole spin on debris 
circularization 

First, we examine the circularization of stellar debris in the 
radiatively efficient cooling case for Models 4-6, in which 
the vector of black hole spin is aligned with z-axis. We note, 
from panel (b2) of Figure]^ that the energy dissipation rate 
of Model 6 is highest among these three models. This shows 
that the enhanced apsidal precession of retrograde black hole 
spin makes debris circularization more efficient than in the 
non-spinning black hole case, whereas prograde black hole 
spin substantially delays debris circularization by decreas¬ 
ing the apsidal precession rate (Merritt et al. ]|2010[ ). Debris 
circularization progresses rapidly but saturates at a specific 
energy slightly less than et in Models 4 and 6 without reach¬ 
ing Cc, at which point debris circularization would be com¬ 
pleted. This incomplete circularization is because relative 
velocities between debris head and tail at self-intersections 
(near apocentre) are too low to produce significant shock 
dissipation. This could originate from the low resolution of 
the simulation (see section 3.4). 

We also simulate coplanar circularization processes, in 
Models 4-5, for a tidally disrupted star around a spin¬ 
ning SMBH in the radiatively inefficient cooling limit. From 
panel (d2), we see that the specific energy increases with 
time after t = 10. As in Models 1-3, this is because the most 
tightly bound gas is preferentially accreted by the black hole, 
increasing the average specific energy of remaining disk mat¬ 
ter. The circularized accretion torus extends to much larger 
radii than in the radiatively efficient case, as some of the 


gas spreads to more loosely bound circular orbits due to 
hydrodynamic effects. Specific angular momentum also in¬ 
creases with time after t = 10, again because of preferential 
accretion. 

The main difference between the different cooling 
regimes is clear from comparing Figure and Figure 
which show a sequence of snapshots for Model 5 in the ra¬ 
diatively efficient and inefficient cooling cases, respectively. 
These figures have the same formats as Figure with each 
panel in chronological order, from panel (a) to (1). Increased 
thermal energy from shock heating causes the debris tem¬ 
perature to rise significantly at early times in the radiatively 
inefficient case, which produces a geometrically thick accre¬ 
tion disk at late times. This is in sharp contrast to the geo¬ 
metrically thin gas structure produced in the radiatively ef¬ 
ficient case of Figure We have also performed radiatively 
inefficient simulations of Models 7-9. These differ from their 
radiatively efficient counterparts in several qualitative ways, 
which we discuss next. 


The effect of the inclination angle between the black 
hole spin vector and the disk angular momentum vector on 
debris circularization can be seen in panel (b3). When the 
black hole rotates retrogradely as in Models 7 and 8, the 
circularization timescale is about ~ 30P*. This is longer 
than the circularization timescale ~ 20P* of the coplanar 
case (Model 6) as can be seen from the red line in panel 
(b2). However, the specific binding energy saturates at ~ Cc, 
which is lower value than that of Model 6. Also, Models 7 
and 8 periodically alternate between rapid energy dissipa¬ 
tion followed by energy conservation, for 14 < t < 40. This is 
because the stellar debris undergoes both relativistic apsidal 
and nodal precession ( Stone et al.||2015 ). The stellar debris 
dissipates orbital energy in shocks due to apsidal precession. 
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but since the debris undergoes nodal precession with each 
pericenter passage, streams that would exactly self-intersect 
in the coplanar case miss each other by a small amount. 
This is the origin of the saturation phase of specific energy 
evolution. Figure shows a sequence of snapshots of the 
tidal disruption process in Model 7 in the radiatively effi¬ 
cient limit. Each panel shows surface densities projected on 
the x-y plane (left panel) and on the y-z plane (right panel) 
for 0 ^ t ^ 40. After tidal disruption, the debris in this run 
circularizes into a thin ring-like structure that slowly rotates 
around the x-axis due to Lense-Thirring precession. As seen 
in the panels (c)-(e), the multiple ring-like structures formed 
by the cumulative effect of nodal precession enable intersec¬ 


tions with large relative velocity even at late times (Stone 
et al.|20T5 ). This enhances both the energy dissipation rate 


and the total amount of dissipated energy above those of 
the coplanar case. 


4.2 Nodal precession caused by the 
Lense-Thirring effect 

Figure shows the evolution of a tilt angle ^tiit and preces¬ 
sion angle Oprec (equivalent to the nodal angle in standard 
orbital elements) for Models 7-9. These angles are defined 
as ( Nelson Papaloizou|2000l Fragile Anninos|2005 ) 


^tiit = arccos 


Or. 


= arccos 


Jbh ' Jd 

|«/bh • Jd\ 

«/bH X Jd 
|JbH X Jd| 


y 


(19) 


( 20 ) 


where Jd, Jbh and y are the angular momentum vector of 
the stellar debris, black hole spin vector, and the unit vector 
of the y-axis, respectively. Here, we set Jbh = (1,0,0) for 
Model 7, Jbh = (l/v^, 0, l/\/2) for Model 8, and Jbh = 
(-1/V2,0,-1/V2) for Model 9. 

We see from this figure that the tilt angle remains con¬ 
stant, but the precession angle increases with time in all 
three models. This is in rough agreement with our expecta¬ 
tions for nodal precession caused by frame dragging torques. 
Let us focus on the radiatively efficient case of Model 7 as 
a concrete example. In the early stages of debris circular¬ 
ization, the precession angle per orbital period in our simu¬ 
lations is roughly ~ O.OItt. This is in good agreement with 
the theoretically expected (to lowest PN order) change per 
orbit of the nodal angle for a test particle orbiting around a 
spinning black hole. 


Obt,^ 

27T 






?'p(l + e,) 


3/2 


V2\rsJ 


-3/2 


Vl + eJ 


3/2 


,(21) 


where Vp = a*(l — e*) is the pericentre distance. This is 


given by equations (6b) and (8) of Merritt et al. (2010), 


using equation @ of this paper, which shows the changes 
per orbit of the nodal angle due to frame dragging torques. 

In the late stages of debris circularization, stellar debris 
has begun to form an accretion disk. In this phase, the pre¬ 
cession angle per orbit grows to ~ 0.027r. This is in good 
agreement with the the change per orbit of the nodal angle 
of an accretion disk, which is given by 


^LT,disk _ 

27r v/2 Vr-s 


r'’m 


2C-5\ 1 - 


2C + 1 y 1 - (ri/ro)^“^/2 ’ 


( 22 ) 


where n, Tq, C Ihe inner radius, outer radius, and power 
law index of the disk surface density, respectively. We have 
approximated (^LT,disk/27r) = P*/tLT,disk, where t LT,disk is 
the precession timescale given by equation (3) of Stone & 
Loeb|(|2012a|). 


The local Lense-Thirring precession timescale is given 


by 


tlt 


^=2- - 




(23) 


where Qlt = 2x((JM)^/(rc)^ ( Bardeen Petterson||l975 ). 
The vertical viscosity timescale is 


2^2 


1 


miTiv) 


(24) 


where y = 2(l + 7ass)/(<^ss(4 + CKss)) is the ratio of the ver- 


tical viscosity to horizontal viscosity _ 

& Price] 2010). For ass ^ 1, ?7 ~ l/2ags ([Papaloizou 


Pringle |1983 ). If the local precession timescale is shorter 


than the vertical viscous timescale, the disc is not warped 
and precess as a rigid body. This condition is given by 

2/3 


IL>(^V (K 


rs 


Stt J 


(?) 


- 4/3 


2/3 2/3 


X 


ss 


(25) 


for ass ^ 1- For typical parameters of the geometrically 
thick disc: H/r ~ 1, ass = 0.1, and y = 0.9, the condi¬ 
tion gives r > 6 X 10“^rs. We obtain the condition r > n 
for the geometrically thin disc case H/r ~ 0.01. While the 
disc rigidly precesses in the radiatively inefficient cases, it is 
unlikely to in the efficient cooling cases at the current low 
resolution simulation. 

The roughly rigid body precession seen in our efficient 
cooling simulations is likely an artifact of the narrow ring¬ 
like configuration that is preserved here because of the short 
simulation runtime. If these simulations were run for longer 
than a viscous time, the spreading gas rings would likely de¬ 


velop Bardeen-Petterson warps (Bardeen & Petterson 1975). 
More realistically, however, the disc viscosity should be 
much less than that of our simulations, because the arti¬ 
ficial viscosity used here would be overestimated by the low 
resolution simulations. Therefore, the disc would be warmer 
and geometrically thicker, causing the disc to be closer to the 
wavelike regime of warp propagation. Even for large, ~ 10® 
solar mass SMBHs where the gas streams and an eccentric 
disc may initially be radiatively efficient (as we argue), a 
fully circularized disc will likely be radiatively inefficient. In 
most TDEs, therefore, the disc can rigidly precess, because 
the Lense-Thirring timescale is longer than the timescale 


the warps propagate as a wave. Franchini et al. (2016) have 


recently shown that the internal dissipation in the disc can 
quench the rigid-body precession while it remains geomet¬ 
rically thick, although typically after the accumulation of 
several precession periods. Our simulations do not run long 
enough to observe such an effect. 

There is a remarkable difference between Models 7-9 in 
the radiatively efficient case, and an even starker difference 
between radiatively efficient and inefficient simulations of 
these models. In Figurethe dashed lines show the evolu¬ 
tion of nodal angles for the geometrically thick disks of our 
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Figure 14. A sequence of snapshots of the tidal disruption process in Model 5 (a* = 5/3, e* = 0.7, /3 = 2, x = 0 - 9 : and z = 0°) in 
the radiatively efficient regime. They run from panel (a) to panel (1) in chronological order. Each panel shows surface densities projected 
on the x-y plane over two orders of magnitude in a logarithmic scale for 0 ^ t ^ 40 (see also color bar), where t is in units of P*. The 
black hole is set at the origin. The run time is annotated at the top-right corner, while the number of SPH particles are indicated at the 
bottom-right corner. The dashed circle indicates the tidal disruption radius. 


radiatively inefficient scenarios, while solid lines show this 
in the radiatively efficient ones. In the radiatively efficient 
case, there are two stages in the nodal angle evolution. In 
the early stage until t ~ 12P*, all the nodal angle curves 
are overlapping as theoretically expected. In the later stage, 
they split in three. These differences between Models 7-9 
are explained by their somewhat different values of n and 
To, which are caused by the difference of formation speed of 
the debris rings. As shown in panel (c2) of Figure]^ the de¬ 
bris ring of Model 8 forms more quickly than that of Model 
7, which forms more quickly than that of Model 9. 

The green solid line of Figure 17 denotes the evolution 
of the nodal angle of the radiatively efficient A^sph = 500K 
simulation for Model 7. The comparison between the lOOK 
and 500K cases shows that the precession rate of the 500K 
case is smaller than the lOOK case in the latter evolution 
stage, because the inner edge radius of the disc is larger 
than the lOOK case by the smaller artificial viscosity owing 
to the smaller smoothing length of the 500K case (i.e. higher 


resolution). It is seen from equation (22) that the larger 
inner edge radius makes the nodal precession slower. 

On the other hand, there is only one stage in the nodal 
angle evolution in the radiatively inefficient limit. Since the 
geometrically thick disks of Models 7-9 have almost identi¬ 
cal surface density profiles in the late stage, all the curves 
overlap for 0 ^ t ^ 40. If we adopt representative values for 
these realistic disks, e.g. Vi = 1.5n, Vo = 4.5n, and C = 2.0 
from Figure the precession timescale is estimated to be 
tLT,disk ~ 2 X 10^ P*. This implies that the nodal angles 
of radiatively inefficient disks increase slowly, and at rates 
decreasing with time. 

After debris circularization, the accretion disk will vis¬ 
cously evolve and accrete onto the black hole. In the efficient 
cooling case, tvis is of the order of 10 ^ P*, as H/r ~ 0 . 01 . 
Since the viscous timescale is clearly longer than the preces¬ 
sion timescale, < lO^P*, the disk will rapidly precess around 
the black hole spin axis as the material accretes onto the 
black hole; this could potentially be imprinted on the ob¬ 
served light curve. In contrast, tvis is of the order of P* in 
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Figure 15. A sequence of snapshots of the tidal disruption process in Model 5 (a* =5/3, e* = 0.7, /3 = 2, x = O-Q, and z = 0°) in the 
radiatively inefficient limit. They run from panel (a) to panel (1) in chronological order. The figure format is the same as Figure [m] 


the inefficient cooling case, where i//r ~ 1 is adopted. The 
viscous timescale is therefore much shorter than the pre¬ 
cession timescale ~ lO^P*. In this regime, which is much 
more realistic for circularized disks in both eccentric and 
parabolic TDEs, most of the disk mass will drain onto the 
SMBH before it can significantly precess. In eccentric TDEs, 
this makes it unlikely that evidence for spin-induced preces¬ 
sion can be found in the observed light curve. In parabolic 
TDEs, the disk is continuously replenished with new mat¬ 
ter, so that unlike the eccentric TDE regime, disk precession 
may be observable. 


5 SUMMARY & DISCUSSION 

We have performed numerical simulations of circularization 
and subsequent accretion disk formation during the tidal dis¬ 
ruption of stars on bound orbits around spinning and non¬ 
spinning SMBHs. We have approximated relativistic effects 
with simple Post-Newtonian corrections up to 2PN, includ¬ 
ing the lowest order spin terms (1.5PN). We have considered 
relatively low orbital eccentricity (e = 0.7 — 0.9) with mod¬ 
est penetration factors /3 = 1 — 2. We have found that debris 


circularization depends crucially on the efficiency of the ra¬ 
diative cooling, and therefore simulated our nine models in 
two limiting regimes: radiatively efficient, and radiatively 
inefficient. Our main conclusions are as follows: 

(i) There are two stages in debris circularization if ra¬ 
diative cooling is inefficient: In the early stage, the stellar 
debris, stretched by tidal disruption, orbits and forms a ge¬ 
ometrically thick ring-like structure around the black hole 
due to shock heating caused by orbital self-intersections in¬ 
duced by relativistic precession. In the late stage, the ring¬ 
like structure rapidly spreads via viscous diffusion, forming 
a geometrically thick accretion disk. In contrast, if radia¬ 
tive cooling is efficient, the stellar debris circularizes into a 
geometrically thin ring-like structure. 

(ii) In relatively low eccentricity tidal disruptions, the 
stellar debris is clearly optically thick and its photon diffu¬ 
sion timescale is longer than the shock heating timescale, so 
stellar debris will circularize as in our radiatively inefficient 
cooling simulations. However, in a parabolic tidal disrup¬ 
tion for Mbh ^ 2 X 1O®M0, our radiatively efficient simu¬ 
lations may be more relevant, because the photon diffusion 
timescale can be shorter than the fallback timescale. 

(iii) In the radiatively efficient regime, debris circularizes 
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Figure 16. A sequence of snapshots of the tidal disruption process in Model 7 (a* = 5/3, 6=,= = 0.7, (3 = 2^ x = —0.9, and i = 90°) in 
the radiatively efficient case. Each panel shows surface densities projected on x-y plane (left panel) and on y-z plane (right panel) for 
0 ^ t ^ 40. The other figure formats are the same as Figure 


more quickly for retrograde spins than for no spin, and more 
quickly for no spin than for prograde spins. This is because 
retrograde spin increases the apsidal shift per orbit, while 
prograde spin decreases it. Increased apsidal precession both 
reduces the time it takes for the debris head to catch its tail 
(increasing circularization rates for eccentric TDEs) and in¬ 
creases the relative velocity at the stream self-intersection 
point (increasing circularization rates for both eccentric and 
parabolic TDEs). This spin dependence is largely absent 
from the radiatively inefficient regime, where the increased 
stream thickness due to heating dominates more subtle GR 
effects. As discussed in section 3.4, our radiatively ineffi¬ 
cient simulations appear numerically well-converged during 
the period of debris circularization, but the radiatively effi¬ 
cient simulations are underestimating the efficiency of shock 
dissipation in lOOK particles’ simulations. 

(iv) When the the black hole spin axis is initially mis¬ 
aligned with the debris angular momentum vector, the cir¬ 
cularized debris ring or disk processes due to Lense-Thirring 
torque. While the tilt angles remains constant during de¬ 
bris circularization in both the radiatively inefficient and 
efficient cooling cases, the nodal precession angles vary with 


time, indicating approximately solid-body precession. While 
nodal angles in both cases evolve in the same manner prior 
to significant energy dissipation in shocks, they evolve much 
more slowly in the radiatively inefficient regime during sub¬ 
sequent dissipative circularization, because the debris forms 
a geometrically thick accretion disk with much larger radial 
extent. 

(v) In the radiatively efficient regime, debris circulariza¬ 
tion is significantly impeded by misaligned SMBH spin. De¬ 
spite this retarding effect, however, it increases the total 
energy dissipated during debris circularization by the end of 
our simulations, although it is unlikely that the radiatively 
efficient cooling regime applies at late times, once debris 
is mostly circularized. Nodal precession of cool, thin debris 
streams may result in significantly reduced shock dissipa¬ 
tion at self-intersection points, as the streams can miss each 
other completely, or suffer only a grazing collision. 

Expanding on conclusion (v), we discuss how applica¬ 
ble these simulations of eccentric TDEs are to the presum¬ 
ably more common parabolic TDE scenario. The dynamics 
of stream self-intersections are very similar: as implied by 
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Figure 17. Evolution of the tilt and nodal angles for Models 7-9 in both the radiatively efficient and inefficient cooling cases. The 
tilt angle ^tilt shows the angle between the debris angular momentum vector and the black hole spin vector, and the nodal angle 6prec 
shows the angle between the debris angular momentum vector and the z-axis. These two angles are averaged over all SPH particles 
and normalized by 27r. The dashed lines show the nodal angles in the radiatively inefficient cooling cases, and the solid lines denote the 
nodal angles in the radiatively efficient cooling cases. The dotted lines denote the tilt angles in both cases. Note that the the green solid 
line denotes the nodal angle in the radiatively efficient A^sph = 500K particles’ simulation for Model 7. The run time t is in units of 
P. = 2w^rflGM ~ 2.8 hr. 


equation (21), the nodal precession rate per orbit depends 
almost entirely on the pericentre distance. Orbit-averaged 
apsidal precession behaves similarly. Debris streams from 
parabolic TDEs will therefore self-intersect at points very 
similar to the self-intersection points in the radiatively effi¬ 
cient simulation of Model 1. Also, in our radiatively efficient 
Models 1-3, we see increasing circularization efficiency with 
increasing eccentricity; this is largely a function of stream 
relative velocity at the self-intersection radius. The circular¬ 
ization efficiency can become very low when this radius is 
comparable to the stream apocentre. 


The greater uncertainty in extrapolating our results to 
the parabolic limit is the structure of the streams, which de¬ 
pends sensitively on internal self-gravity and cooling physics. 
Modeling stream structure in parabolic TDEs goes beyond 
the scope of this paper. However, our two extreme cooling 
regimes have outlined the relevant parameter space: if the 
stream structure is comparable to or thinner than that in 
our constant-entropy simulations, then debris circularization 
will depend sensitively on SMBH spin. If stream structure 
is closer to our radiatively inefficient limit, then SMBH spin 
will have very little effect on debris circularization. 

As a final topic of discussion, we consider where the 
by debris circularization goes. Recently, 
showed that vertical advection of radia¬ 
tion caused by magnetic buoyancy transports energy about 
50 times faster than does photon diffusion in the super- 
Eddington, radiation-pressure dominated accretion regimes 
that characterize realistic TDE disks. The advective cooling 
timescale in the vertical direction, tadvz, is thus estimated to 
be of the order of 10^ s for the typical parameters of equa¬ 
tion Therefore, tcooi = min(tdiff, tadvz) is the cooling 
timescale. 


energy dissipated 


Jiang et al. (2014 


If we optimistically assume tcooi ^ tc, then the energy 
dissipated by shocks during debris circularization will radi¬ 


ate away as thermal emission via optically thick radiative 
cooling. The resultant Lc = Jcmax/tc gives an upper limit 
on the circularization luminosity of parabolic TDEs: 


Lc 

Ledd 


X 



/ Mbh a 
\10^Mq ) 


- 5/2 


(26) 


where we adopt tc = iVperitfb with the number of pericen¬ 
tre passages A^peri during debris circularization, and L^dd = 
47rGMBH^pc/(7T is the Eddington luminosity with the pro¬ 
ton mass rup and the Thomson cross section ctt- In the cases 
of eccentric TDEs or radiatively inefficient parabolic ones, 
this luminosity is reduced by a factor tc/^diff • The actual ra¬ 
diated energy could be further reduced by a transition to a 
radiatively inefficient cooling flow or by a super-Eddington 
outflow ( Strubbe Quataert|20lf ). However, it is clear that 
the circularization energy budget, especially for parabolic 
TDEs, is substantial. Such emission could be observed as a 
precursor to the parabolic TDE flare powered by accretion 
onto the SMBH for Mbh ^ 2 x 10® M©, where the radiatively 
efficient scenario is applicable. For typical parameters, the 
upper limit of the precursor luminosity can be comparable 
to the Eddington luminosity. 


A similar idea was proposed by equation (8) of Bog- 


danovic et al. (2014) in the context of tidal disruption of a 


red giant star on a parabolic orbit. However, it is important 
to evaluate whether the energy dissipated by shocks can ac¬ 
tually be radiated away within a fallback time. If this is not 
the case, then any precursor signal from circularization lu¬ 
minosity will be strongly suppressed as with our radiatively 
inefficient cooling simulations. We have found that the dissi¬ 
pated energy will be efficiently radiated away for black holes 
above a critical mass: ~ 2 x 10®M© for electron scattering 
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Figure 18. A sequence of snapshots of the tidal disruption process in Model 7 (a* = 5/3, 6=,= = 0.7, (3 = 2^ x = —0.9, and i = 90°) in 
the radiatively inefficient limit. Each panel shows surface densities projected on x-y plane (left panel) and on y-z plane (right panel) for 
0 ^ t ^ 40. The other figure formats are the same as Figure 


opacity. The precursor luminosity can be more accurately 
obtained by self-consistently calculating the cooling rate of 
the emitted region of the debris and its stream structure. A 
more detailed study will be done in the future. 


In summary, eccentric TDEs serve as a valuable and 
computationally tractable testbed for the physics of circu¬ 
larization, which is extremely challenging to simulate for 
parabolic TDEs around SMBHs. Our first ever simulations 
of stellar tidal disruption around spinning SMBHs confirm 
past analytic predictions of disk precession, and indicate 
that black hole spin may imprint itself, via circularization 
delays, onto mass fallback rates in parabolic TDEs. Future 
study of debris stream dynamics in the parabolic limit is 
required to make more definite predictions, but for now we 
conclude that SMBH spin is of crucial importance for the 
circularization of thin debris streams, a process aided by ret¬ 
rograde or aligned spin, and hindered by prograde or mis¬ 
aligned spin. 


APPENDIX A: SIMPLE TREATMENT OF 
RELATIVISTIC EFFECTS IN SPH 

Here we briefly describe how to incorporate PN correction 
terms relevant for the tidal disruption problem into an exist¬ 
ing SPH code. Specifically, we add the PN acceleration terms 
into the SPH momentum equation, and leave the other SPH 
equations unchanged. 

The SPH momentum equation for the z-th particle with 
PN corrections can be written by 


+ Yiij I Vj VP (rij , hij ) 


H -1.5PN H -jCti 2PN, (A1) 

where Pj and pj are the pressure and density of particle j, 
respectively, and Uij is the standard form of the artificial 


dvi 

dt 


(Pi Pj 

Y Pi 

GM(rij) Vij 


- E 


+ OPN H- TjCii,IPN 
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viscosity ( Bate||1995 ): 


{ — aspuCsinj + Pspunlj)/Pij '^ij ' 'f^ij ^ 0 


0 


• ra > 0 . 


(A 2 ) 


2 ^ \ ^ 2 
- -(Vi-VBH) + ^VBH 


TliBH 


Here, asPH and /3sph are linear and nonlinear artificial vis¬ 
cosity parameters, respectively; pij = (pi + pAI^, vij — 
Vi - Vj, and pij = hijVij • rij^nj + pij) with pfj = O.Olh^j. 
Note that we adopt the standard values, asPH = 1 and 
/^SPH = 2 , for the artificial viscosity parameters in our sim¬ 
ulations. The artificial viscosity consists of two terms: the 
first term that is linear in the velocity differences between 
particles, which produces a shear and bulk viscosity, and 
the second term that is quadratic in the velocity differences, 
which is needed to eliminate particle interpenetration in high 
Mach number shocks. W is the weighting function called by 
kernel. We adopt the standard cubic-spline kernel for 3D 
( Bide|[T^ : 


GM^ 
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if 1 ^ s < 2 (A3) 

otherwise, 
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—^(n^BH • Vi) 


+ — (n^BH'^BH) 
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The second term of right hand side of equation (Al) 


shows the self-gravitating force acting on particle i from all 
the other SPH particles, where 


M{nj) = 


nr 

j) = 4n 

Jo 


p{r) dr. 


The third term to the sixth term show the gravitational 
forces acting on particle i from the black hole. The Newto¬ 
nian (OPN), IPN, and 2PN accelerations among them can 
be written by ( Blanchet||2QQ6 ) 


ai,oPN = — - 


GMbb 
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ViBH(A4) 


where vbb and vbh are the position and velocity vector of 
the black hole particle, respectively, and t^bh = Vi — vbb, 
ViBB = Vi — vbh, and riiBB = {vi — 'rBH)/|T’z — rBH|. The 
1.5PN acceleration due to the black hole spin can be also 
written by Faye et al. (2007) as 


tti,1.5PN 


GMbb 
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Si • (n^BH X VzBh) 
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Sbb • (n^BH X Vzbh) 


Mbh 

+ 3(n zBH • '^iBH) 
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Mbb 

ViBB X Sbb 
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(A5) 


where Si and Sbb are the spin vectors of particle i and 
black hole, respectively. In our smulation. Si = 0 and 
*S'bh = GM^ylX^ are adopted, where y is the black hole 
spin parameter with value of 0 ^ y ^ 1 and s is the unit 
vector of spin angular momentum of the black hole. 


Al Post-Newtonian corrections to binding energy 
and angular momentum 

The binding energy and angular momentum of SPH particle- 
sare also corrected by PN approximations. The total binding 
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energy of the system is given by 

- ^SPH 

£=^E,, 

i=l 


(A 6 ) 
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where the binding energy of a SPH particle, Ei, can be cor¬ 
rected ( Blanchet]|2006 Faye et aT]|2Q07 ) as 


Ei — F^z, 0PN H-^^i,lPN H-^^i,1.5PN H-j^i,2PN, 


(A7) 


where 

^i,OPN 

^i,lPN 




1 / 2 , ,. 2 N GrriiMBB 

-{rriiVi + Mbh^bh)-, 

2 TiBB 

G^tti^Mbb rriivf ^ GuiiMBB 


zBH 


^BH 


1 3 

— j(niBH • Vi){niBB • Vbb) + -Vi 


- -{Vi-VBB) 


G^MBB^^i Mbh^bh GMBurni 


iBH 


^bh 


1 3 

— -(n^BH • VBH)(niBH • Vi) + 2 '^BH 


- -rivBB'Vi) 


i,1.5PN = 


^i,2PN = 


4 

GMbb 


[Si • (riiBH X Vi)] 


Girii r ^ / N] 

—[Sbh • (niBH X vbh)J, 


G^mfMBH 5 e 

- + TTlVriiVi 


2 ^?BH 




16 


G^ MBB^^i \ ^ i\/r 6 

- a:- * 


' zBH 

G^ttiIMbb 


29. .2 

^(riiBH • Vi) 


13 


(riiBH • Vi)(niBH • vbh) 


^(^zBH • '^Bh)^ - -vl^ + J^BH 


AfgjjTTli 
' iBH 


+ —{riiBB • Vbb)^ (riiBB • Vi)‘^ 

— -(^iBH • 'l^BH)(TliBH ' 'l^O'^BH 

o 

13, .2 2 21 4 

“ -^(^zBH • Vi) Vbh + “^^BH 

13 

+ — (riiBH • 'yBH)^('yz • Vbh) 


+ j(niBH • 'yBH)(niBH • Vi)(Vi • vbh) 

~ • '^Bh) + — • Vbh) E Y^ViVBB 

The total angular momentum vector of the system can 
be written by 


-^SPH 

Z=1 


(A 8 ) 


where the angular momentum vector of a SPH particle, Ji, 


is corrected by de Andrade et ah (2001); Faye et ah (2007) 
as 

Ji — «/z,0PN + ^«/z,lPN + ^«/z,1.5PN + ^«/z,2PN. (^9) 

Here, each term of the right-hand side can be written as the 
ingredient-label format by 

'7i,0PN = Slmn{vflir^M bbTbbVbb) 


J 


z,lPN — ^Imn 


m n I SGrUiMBB , rUiV^ 

Eg Vi -+ 

riBH 2 


m n ’^GlTlilElBB rn n GlTlilElBB / \ 

Vi ^BH—^r::- \-ri -(riiBH-Vi) 


2riBH 


,■ 


^z,2PN — 


J m n TGlTlilElBB / \ 

^lmn\ Vi Vbh ~a ('^zBH " Vi) 


bG^mtMBB ^ IG^tthMIb 


+ 


29. .2 13. .. . 

— V^iBH-'l^BH)-^ V^^BH • 'l^BH)(niBH • 

1 / n 2 3 2 7 2 , GrriiMBB 

-(riiBB'Vi) — -vbbE -Vi H- 

2 Z 4 riBH 


-(riiBH • Vi)^{niBB • Vbh) 


+ — (riiBH • (riiBH • 'I^Bh)^ 

16 

9/ \ / \2l3/ \22 

— -(riiBH • r^i)(riiBH • r^Bnjr’i-^(riiBH • i^bh) Vi 

o o 

21 4 13, ,2/ X 

+ E — (riiBH • Vi) (Vi • Vbh) 

o o 

3 

+ y (riiBH • 'yz)(^^BH • VBH)('yz • I^Bh) 


4^-?bh 


riBH 


0'i^2 

zBH 


— ^(^^bh • i^bh)^ + -Vi 


— 4(vi • vbh) + 2ve 


I m n 

E Vi Vbh 


7 G^ TTli Mbi 

4 T 2 

zBH 


GrriiMBB / 1 / x2 

- -y (riiBH • Vi) 

riBH \ 8 

1 13 

y(niBH • Vz)(niBH • VBh) + —(riiBH • vbh) 

4 o 


+ ^C’i- ■^Bh) - y^’BH 


, rn Ti 

E Vi Tbh 


G rrii Mbb 


' zBH 


I 29. \ I \ \ I GrriiMBB 

X ( — — (riiBH • Vi) + - (riiBH • i^bh) I H- 2 - 

4 4 / r -T^TT 


X ( — -(niBH • Vi)^ — -(riiBH • Vi)^{niBH ■ vbh) 
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+ -{riiBii ■ Vi)vi -{niBu ■ VBii)vi 


— -^{riiBB ■ Vi){Vi ■ Vbu) 


(AlO) 


where eimn shows Levi-Civita tensor with {l,m,n = 1, 2, 3), 
and the spin term can be written as a vector format by 


«/i,1.5PN = 


+ 

+ 

+ 


GMbh / c* ^ 

-2-^ {'TliBH X Si) 

^iBH 

X (riiBB X Sbk) 

^zBH 

GMbh r^/ \ o i i ^ 2 c* 

-[2(niBH • Si) • riiBH — + -Vi Si 

TiBH ^ 

(vi • Si)vi 

-^[2(niBH • *S'bh) • 'H-iBH — *S'bh] + -^BH'^BH 

TiBH 2 


(vbh • »S'bh)'I^BH. 


(All) 
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